Load the necessary libraries
library(car) # for regression diagnostics
library(broom) # for tidy output
library(ggfortify) # for model diagnostics
library(sjPlot) # for outputs
library(knitr) # for kable
library(effects) # for partial effects plots
library(emmeans) # for estimating marginal means
library(MASS) # for glm.nb
library(MuMIn) # for AICc
library(tidyverse) # for data wrangling
library(brms)
library(tidybayes)
library(broom.mixed) # for tidying MCMC outputs
library(patchwork) # for multiple plots
library(standist) # for visualizing distributions
library(rstanarm)
library(ggeffects)
library(bayesplot)
library(rstan)
library(DHARMa)
source("helperFunctions.R")
Starlings
Format of starling_full.RSV data files
| SITUATION | MONTH | MASS | BIRD |
|---|---|---|---|
| tree | Nov | 78 | tree1 |
| .. | .. | .. | .. |
| nest-box | Nov | 78 | nest-box1 |
| .. | .. | .. | .. |
| inside | Nov | 79 | inside1 |
| .. | .. | .. | .. |
| other | Nov | 77 | other1 |
| .. | .. | .. | .. |
| tree | Jan | 85 | tree1 |
| .. | .. | .. | .. |
| SITUATION | Categorical listing of roosting situations (tree, nest-box, inside or other) |
| MONTH | Categorical listing of the month of sampling. |
| MASS | Mass (g) of starlings. |
| BIRD | Categorical listing of individual bird repeatedly sampled. |
starling <- read_csv("../public/data/starling_full.csv", trim_ws = TRUE)
starling %>% glimpse()
## Rows: 80
## Columns: 5
## $ MONTH <chr> "Nov", "Nov", "Nov", "Nov", "Nov", "Nov", "Nov", "Nov", "No…
## $ SITUATION <chr> "tree", "tree", "tree", "tree", "tree", "tree", "tree", "tr…
## $ subjectnum <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1…
## $ BIRD <chr> "tree1", "tree2", "tree3", "tree4", "tree5", "tree6", "tree…
## $ MASS <dbl> 78, 88, 87, 88, 83, 82, 81, 80, 80, 89, 78, 78, 85, 81, 78,…
Model formula: \[ y_i \sim{} \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i = \beta_0 + \boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i}\\ \boldsymbol{\gamma} = \gamma_0\\ \beta_0 \sim{} \mathcal{N}(0, 100)\\ \beta \sim{} \mathcal{N}(0, 10)\\ \gamma_0 \sim{} \mathcal{N}(0, \sigma_1^2)\\ \sigma \sim{} \mathcal{cauchy}(0, 2)\\ \sigma_1 \sim{} \mathcal{cauchy}(0, 2)\\ \]
where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of roosting situation and month on starling mass. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual birds.
starling <- starling %>%
mutate(
BIRD = factor(BIRD),
SITUATION = factor(SITUATION),
MONTH = factor(MONTH, levels = c("Nov", "Jan"))
)
ggplot(starling, aes(y = MASS, x = MONTH)) +
geom_boxplot() +
facet_grid(~SITUATION)
## Better still
ggplot(starling, aes(y = MASS, x = MONTH, group = BIRD)) +
geom_point() +
geom_line() +
facet_grid(~SITUATION)
Conclusions:
In rstanarm, the default priors are designed to be
weakly informative. They are chosen to provide moderate regularisation
(to help prevent over-fitting) and help stabilise the computations.
starling.rstanarm <- stan_glmer(MASS ~ MONTH * SITUATION + (1 | BIRD),
data = starling,
family = gaussian(),
iter = 5000,
warmup = 2000,
chains = 3,
thin = 5,
refresh = 0
)
starling.rstanarm %>% prior_summary()
## Priors for model '.'
## ------
## Intercept (after predictors centered)
## Specified prior:
## ~ normal(location = 84, scale = 2.5)
## Adjusted prior:
## ~ normal(location = 84, scale = 17)
##
## Coefficients
## Specified prior:
## ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
## Adjusted prior:
## ~ normal(location = [0,0,0,...], scale = [33.25,38.40,38.40,...])
##
## Auxiliary (sigma)
## Specified prior:
## ~ exponential(rate = 1)
## Adjusted prior:
## ~ exponential(rate = 0.15)
##
## Covariance
## ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details
This tells us: - for the intercept (when the family is Gaussian), a
normal prior with a mean of 84 and a standard deviation of 17 is used.
The 2.5 is used for scaling all parameter standard deviations. The value
of 84 is based on the mean of the response
(mean(starling$MASS)) and the scaled standard deviation of
17 is based on multiplying the scaling factor by the standard deviation
of the response.
2.5 * sd(starling$MASS)
## [1] 16.73225
2.5 * sd(starling$MASS) / apply(model.matrix(~ MONTH * SITUATION, starling)[, -1], 2, sd)
## MONTHJan SITUATIONnest-box
## 33.25470 38.39922
## SITUATIONother SITUATIONtree
## 38.39922 38.39922
## MONTHJan:SITUATIONnest-box MONTHJan:SITUATIONother
## 50.27638 50.27638
## MONTHJan:SITUATIONtree
## 50.27638
rstanarm, this is a exponential with a rate
of 1 which is then adjusted by division with the standard deviation of
the response.1 / sd(starling$MASS)
## [1] 0.149412
Lets now run with priors only so that we can explore the range of values they allow in the posteriors.
starling.rstanarm1 <- update(starling.rstanarm, prior_PD = TRUE)
starling.rstanarm1 %>%
ggpredict(~ MONTH * SITUATION) %>%
plot(add.data = TRUE)
Conclusions:
The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]
When defining our own priors, we typically do not want them to be scaled.
If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucked out of thin air):
mean(starling$MASS)2.5*sd(starling$MASS)2.5*sd(starling$MASS) / apply(model.matrix(~MONTH*SITUATION, starling)[,-1], 2, sd)2.5*sd(starling$MASS) / apply(model.matrix(~MONTH*SITUATION, starling)[,-1], 2, sd)2.5*sd(starling$MASS) / apply(model.matrix(~MONTH*SITUATION, starling)[,-1], 2, sd)1 / sd(starling$MASS)I will also overlay the raw data for comparison.
starling.rstanarm2 <- stan_glmer(MASS ~ MONTH * SITUATION + (1 | BIRD),
data = starling,
family = gaussian(),
prior_intercept = normal(84, 17, autoscale = FALSE),
prior = normal(0, c(33, 39, 39, 39, 50, 50, 50), autoscale = FALSE),
prior_aux = rstanarm::exponential(0.15, autoscale = FALSE),
prior_covariance = decov(1, 1, 1, 1),
prior_PD = TRUE,
iter = 5000,
warmup = 1000,
chains = 3,
thin = 5,
refresh = 0
)
starling.rstanarm2 %>%
ggpredict(~ SITUATION * MONTH) %>%
plot(add.data = TRUE)
Now lets refit, conditioning on the data.
starling.rstanarm3 <- update(starling.rstanarm2, prior_PD = FALSE)
posterior_vs_prior(starling.rstanarm3,
color_by = "vs", group_by = TRUE,
facet_args = list(scales = "free_y")
)
Conclusions:
ggemmeans(starling.rstanarm3, ~ SITUATION * MONTH) %>% plot(add.data = TRUE)
ggpredict(starling.rstanarm3, ~ SITUATION * MONTH) %>% plot(add.data = TRUE)
In brms, the default priors are designed to be weakly
informative. They are chosen to provide moderate regularisation (to help
prevent over fitting) and help stabilise the computations.
Unlike rstanarm, brms models must be
compiled before they start sampling. For most models, the compilation of
the stan code takes around 45 seconds.
starling.form <- bf(MASS ~ MONTH*SITUATION+(1|BIRD),
family = gaussian()
)
options(width=150)
starling.form %>% get_prior(data = starling)
## prior class coef group resp dpar nlpar bound source
## (flat) b default
## (flat) b MONTHJan (vectorized)
## (flat) b MONTHJan:SITUATIONnestMbox (vectorized)
## (flat) b MONTHJan:SITUATIONother (vectorized)
## (flat) b MONTHJan:SITUATIONtree (vectorized)
## (flat) b SITUATIONnestMbox (vectorized)
## (flat) b SITUATIONother (vectorized)
## (flat) b SITUATIONtree (vectorized)
## student_t(3, 84, 5.9) Intercept default
## student_t(3, 0, 5.9) sd default
## student_t(3, 0, 5.9) sd BIRD (vectorized)
## student_t(3, 0, 5.9) sd Intercept BIRD (vectorized)
## student_t(3, 0, 5.9) sigma default
options(width=80)
The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]
When defining our own priors, we typically do not want them to be scaled.
If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucked out of thin air):
Note, for hierarchical models, the model will tend to want to have a large \(sigma\) in order to fit the data better. It is a good idea to regularise this tendency by applying a prior that has most mass around zero. Suitable candidates include:
I will also overlay the raw data for comparison.
starling %>%
group_by(SITUATION, MONTH) %>%
summarise(
median(MASS),
mad(MASS)
)
standist::visualize("normal(84,20)", xlim = c(0, 200))
standist::visualize("student_t(3, 0, 5.9)",
"gamma(2,0.5)",
"cauchy(0,1)",
xlim = c(-10, 25)
)
priors <- prior(normal(84, 5), class = "Intercept") +
prior(normal(0, 5), class = "b") +
prior(gamma(2, 0.5), class = "sigma") +
prior(cauchy(0, 1), class = "sd")
starling.form <- bf(MASS ~ MONTH * SITUATION + (1 | BIRD),
family = gaussian()
)
starling.brm2 <- brm(starling.form,
data = starling,
prior = priors,
sample_prior = "yes",
iter = 5000,
warmup = 1000,
chains = 3,
thin = 5,
refresh = 0
)
starling.form <- bf(MASS ~ MONTH * SITUATION + (MONTH | BIRD),
family = gaussian()
)
starling.brm3 <- brm(starling.form,
data = starling,
prior = priors,
sample_prior = "yes",
iter = 5000,
warmup = 1000,
chains = 3,
thin = 5,
refresh = 0,
control = list(adapt_delta = 0.99)
)
(l.1 <- starling.brm2 %>% loo())
##
## Computed from 2400 by 80 log-likelihood matrix
##
## Estimate SE
## elpd_loo -233.0 5.1
## p_loo 10.8 1.3
## looic 466.1 10.2
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 74 92.5% 1088
## (0.5, 0.7] (ok) 6 7.5% 642
## (0.7, 1] (bad) 0 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
##
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.
(l.2 <- starling.brm3 %>% loo())
##
## Computed from 2400 by 80 log-likelihood matrix
##
## Estimate SE
## elpd_loo -232.6 5.1
## p_loo 15.4 1.8
## looic 465.2 10.3
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 54 67.5% 825
## (0.5, 0.7] (ok) 22 27.5% 196
## (0.7, 1] (bad) 4 5.0% 1021
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
loo_compare(l.1, l.2)
## elpd_diff se_diff
## . 0.0 0.0
## . -0.4 0.5
starling.brm3 %>% get_variables()
## [1] "b_Intercept" "b_MONTHJan"
## [3] "b_SITUATIONnestMbox" "b_SITUATIONother"
## [5] "b_SITUATIONtree" "b_MONTHJan:SITUATIONnestMbox"
## [7] "b_MONTHJan:SITUATIONother" "b_MONTHJan:SITUATIONtree"
## [9] "sd_BIRD__Intercept" "sd_BIRD__MONTHJan"
## [11] "cor_BIRD__Intercept__MONTHJan" "sigma"
## [13] "r_BIRD[inside1,Intercept]" "r_BIRD[inside10,Intercept]"
## [15] "r_BIRD[inside2,Intercept]" "r_BIRD[inside3,Intercept]"
## [17] "r_BIRD[inside4,Intercept]" "r_BIRD[inside5,Intercept]"
## [19] "r_BIRD[inside6,Intercept]" "r_BIRD[inside7,Intercept]"
## [21] "r_BIRD[inside8,Intercept]" "r_BIRD[inside9,Intercept]"
## [23] "r_BIRD[nest-box1,Intercept]" "r_BIRD[nest-box10,Intercept]"
## [25] "r_BIRD[nest-box2,Intercept]" "r_BIRD[nest-box3,Intercept]"
## [27] "r_BIRD[nest-box4,Intercept]" "r_BIRD[nest-box5,Intercept]"
## [29] "r_BIRD[nest-box6,Intercept]" "r_BIRD[nest-box7,Intercept]"
## [31] "r_BIRD[nest-box8,Intercept]" "r_BIRD[nest-box9,Intercept]"
## [33] "r_BIRD[other1,Intercept]" "r_BIRD[other10,Intercept]"
## [35] "r_BIRD[other2,Intercept]" "r_BIRD[other3,Intercept]"
## [37] "r_BIRD[other4,Intercept]" "r_BIRD[other5,Intercept]"
## [39] "r_BIRD[other6,Intercept]" "r_BIRD[other7,Intercept]"
## [41] "r_BIRD[other8,Intercept]" "r_BIRD[other9,Intercept]"
## [43] "r_BIRD[tree1,Intercept]" "r_BIRD[tree10,Intercept]"
## [45] "r_BIRD[tree2,Intercept]" "r_BIRD[tree3,Intercept]"
## [47] "r_BIRD[tree4,Intercept]" "r_BIRD[tree5,Intercept]"
## [49] "r_BIRD[tree6,Intercept]" "r_BIRD[tree7,Intercept]"
## [51] "r_BIRD[tree8,Intercept]" "r_BIRD[tree9,Intercept]"
## [53] "r_BIRD[inside1,MONTHJan]" "r_BIRD[inside10,MONTHJan]"
## [55] "r_BIRD[inside2,MONTHJan]" "r_BIRD[inside3,MONTHJan]"
## [57] "r_BIRD[inside4,MONTHJan]" "r_BIRD[inside5,MONTHJan]"
## [59] "r_BIRD[inside6,MONTHJan]" "r_BIRD[inside7,MONTHJan]"
## [61] "r_BIRD[inside8,MONTHJan]" "r_BIRD[inside9,MONTHJan]"
## [63] "r_BIRD[nest-box1,MONTHJan]" "r_BIRD[nest-box10,MONTHJan]"
## [65] "r_BIRD[nest-box2,MONTHJan]" "r_BIRD[nest-box3,MONTHJan]"
## [67] "r_BIRD[nest-box4,MONTHJan]" "r_BIRD[nest-box5,MONTHJan]"
## [69] "r_BIRD[nest-box6,MONTHJan]" "r_BIRD[nest-box7,MONTHJan]"
## [71] "r_BIRD[nest-box8,MONTHJan]" "r_BIRD[nest-box9,MONTHJan]"
## [73] "r_BIRD[other1,MONTHJan]" "r_BIRD[other10,MONTHJan]"
## [75] "r_BIRD[other2,MONTHJan]" "r_BIRD[other3,MONTHJan]"
## [77] "r_BIRD[other4,MONTHJan]" "r_BIRD[other5,MONTHJan]"
## [79] "r_BIRD[other6,MONTHJan]" "r_BIRD[other7,MONTHJan]"
## [81] "r_BIRD[other8,MONTHJan]" "r_BIRD[other9,MONTHJan]"
## [83] "r_BIRD[tree1,MONTHJan]" "r_BIRD[tree10,MONTHJan]"
## [85] "r_BIRD[tree2,MONTHJan]" "r_BIRD[tree3,MONTHJan]"
## [87] "r_BIRD[tree4,MONTHJan]" "r_BIRD[tree5,MONTHJan]"
## [89] "r_BIRD[tree6,MONTHJan]" "r_BIRD[tree7,MONTHJan]"
## [91] "r_BIRD[tree8,MONTHJan]" "r_BIRD[tree9,MONTHJan]"
## [93] "prior_Intercept" "prior_b"
## [95] "prior_sigma" "prior_sd_BIRD"
## [97] "prior_cor_BIRD" "lp__"
## [99] "accept_stat__" "stepsize__"
## [101] "treedepth__" "n_leapfrog__"
## [103] "divergent__" "energy__"
starling.brm3 %>%
hypothesis("MONTHJan=0") %>%
plot()
starling.brm3 %>%
posterior_samples() %>%
dplyr::select(-`lp__`) %>%
pivot_longer(everything(), names_to = "key") %>%
filter(!str_detect(key, "^r")) %>%
mutate(
Type = ifelse(str_detect(key, "prior"), "Prior", "Posterior"),
## Class = ifelse(str_detect(key, 'Intercept'), 'Intercept',
## ifelse(str_detect(key, 'b'), 'b', 'sigma')),
Class = case_when(
str_detect(key, "(^b|^prior).*Intercept$") ~ "Intercept",
str_detect(key, "b_SITUATION.*|b_MONTH.*|prior_b") ~ "TREATMENT",
str_detect(key, "sd") ~ "sd",
str_detect(key, "^cor|prior_cor") ~ "cor",
str_detect(key, "sigma") ~ "sigma"
),
Par = str_replace(key, "b_", "")
) %>%
ggplot(aes(x = Type, y = value, color = Par)) +
stat_pointinterval(position = position_dodge(), show.legend = FALSE) +
facet_wrap(~Class, scales = "free")
In addition to the regular model diagnostics checking, for Bayesian analyses, it is also necessary to explore the MCMC sampling diagnostics to be sure that the chains are well mixed and have converged on a stable posterior.
There are a wide variety of tests that range from the big picture, overall chain characteristics to the very specific detailed tests that allow the experienced modeller to drill down to the very fine details of the chain behaviour. Furthermore, there are a multitude of packages and approaches for exploring these diagnostics.
The bayesplot package offers a range of MCMC diagnostics
as well as Posterior Probability Checks (PPC), all of which have a
convenient plot() interface. Lets start with the MCMC
diagnostics.
available_mcmc()
## bayesplot MCMC module:
## mcmc_acf
## mcmc_acf_bar
## mcmc_areas
## mcmc_areas_data
## mcmc_areas_ridges
## mcmc_areas_ridges_data
## mcmc_combo
## mcmc_dens
## mcmc_dens_chains
## mcmc_dens_chains_data
## mcmc_dens_overlay
## mcmc_hex
## mcmc_hist
## mcmc_hist_by_chain
## mcmc_intervals
## mcmc_intervals_data
## mcmc_neff
## mcmc_neff_data
## mcmc_neff_hist
## mcmc_nuts_acceptance
## mcmc_nuts_divergence
## mcmc_nuts_energy
## mcmc_nuts_stepsize
## mcmc_nuts_treedepth
## mcmc_pairs
## mcmc_parcoord
## mcmc_parcoord_data
## mcmc_rank_hist
## mcmc_rank_overlay
## mcmc_recover_hist
## mcmc_recover_intervals
## mcmc_recover_scatter
## mcmc_rhat
## mcmc_rhat_data
## mcmc_rhat_hist
## mcmc_scatter
## mcmc_trace
## mcmc_trace_data
## mcmc_trace_highlight
## mcmc_violin
Of these, we will focus on:
pars <- starling.brm3 %>% get_variables()
pars <- pars %>%
str_extract("^b.Intercept|^b_SITUTATION.*|^b_MONTH.*|[sS]igma|^sd.*") %>%
na.omit()
pars
## [1] "b_Intercept" "b_MONTHJan"
## [3] "b_MONTHJan:SITUATIONnestMbox" "b_MONTHJan:SITUATIONother"
## [5] "b_MONTHJan:SITUATIONtree" "sd_BIRD__Intercept"
## [7] "sd_BIRD__MONTHJan" "sigma"
## [9] "sigma"
## attr(,"na.action")
## [1] 3 4 5 11 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
## [20] 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
## [39] 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
## [58] 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84
## [77] 85 86 87 88 89 90 91 92 93 94 96 97 98 99 100 101 102 103 104
## attr(,"class")
## [1] "omit"
starling.brm3 %>% mcmc_plot(type = "trace", pars = pars)
The chains appear well mixed and very similar
starling.brm3 %>% mcmc_plot(type = "acf_bar", pars = pars)
There is no evidence of auto-correlation in the MCMC samples
starling.brm3 %>% mcmc_plot(type = "rhat_hist")
All Rhat values are below 1.05, suggesting the chains have converged.
neff_hist (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
starling.brm2 %>% mcmc_plot(type = "neff_hist")
Ratios all very high.
starling.brm3 %>% mcmc_plot(type = "combo", pars = pars)
starling.brm3 %>% mcmc_plot(type = "violin", pars = pars)
The rstan package offers a range of MCMC diagnostics.
Lets start with the MCMC diagnostics.
Of these, we will focus on:
starling.brm3 %>% get_variables()
## [1] "b_Intercept" "b_MONTHJan"
## [3] "b_SITUATIONnestMbox" "b_SITUATIONother"
## [5] "b_SITUATIONtree" "b_MONTHJan:SITUATIONnestMbox"
## [7] "b_MONTHJan:SITUATIONother" "b_MONTHJan:SITUATIONtree"
## [9] "sd_BIRD__Intercept" "sd_BIRD__MONTHJan"
## [11] "cor_BIRD__Intercept__MONTHJan" "sigma"
## [13] "r_BIRD[inside1,Intercept]" "r_BIRD[inside10,Intercept]"
## [15] "r_BIRD[inside2,Intercept]" "r_BIRD[inside3,Intercept]"
## [17] "r_BIRD[inside4,Intercept]" "r_BIRD[inside5,Intercept]"
## [19] "r_BIRD[inside6,Intercept]" "r_BIRD[inside7,Intercept]"
## [21] "r_BIRD[inside8,Intercept]" "r_BIRD[inside9,Intercept]"
## [23] "r_BIRD[nest-box1,Intercept]" "r_BIRD[nest-box10,Intercept]"
## [25] "r_BIRD[nest-box2,Intercept]" "r_BIRD[nest-box3,Intercept]"
## [27] "r_BIRD[nest-box4,Intercept]" "r_BIRD[nest-box5,Intercept]"
## [29] "r_BIRD[nest-box6,Intercept]" "r_BIRD[nest-box7,Intercept]"
## [31] "r_BIRD[nest-box8,Intercept]" "r_BIRD[nest-box9,Intercept]"
## [33] "r_BIRD[other1,Intercept]" "r_BIRD[other10,Intercept]"
## [35] "r_BIRD[other2,Intercept]" "r_BIRD[other3,Intercept]"
## [37] "r_BIRD[other4,Intercept]" "r_BIRD[other5,Intercept]"
## [39] "r_BIRD[other6,Intercept]" "r_BIRD[other7,Intercept]"
## [41] "r_BIRD[other8,Intercept]" "r_BIRD[other9,Intercept]"
## [43] "r_BIRD[tree1,Intercept]" "r_BIRD[tree10,Intercept]"
## [45] "r_BIRD[tree2,Intercept]" "r_BIRD[tree3,Intercept]"
## [47] "r_BIRD[tree4,Intercept]" "r_BIRD[tree5,Intercept]"
## [49] "r_BIRD[tree6,Intercept]" "r_BIRD[tree7,Intercept]"
## [51] "r_BIRD[tree8,Intercept]" "r_BIRD[tree9,Intercept]"
## [53] "r_BIRD[inside1,MONTHJan]" "r_BIRD[inside10,MONTHJan]"
## [55] "r_BIRD[inside2,MONTHJan]" "r_BIRD[inside3,MONTHJan]"
## [57] "r_BIRD[inside4,MONTHJan]" "r_BIRD[inside5,MONTHJan]"
## [59] "r_BIRD[inside6,MONTHJan]" "r_BIRD[inside7,MONTHJan]"
## [61] "r_BIRD[inside8,MONTHJan]" "r_BIRD[inside9,MONTHJan]"
## [63] "r_BIRD[nest-box1,MONTHJan]" "r_BIRD[nest-box10,MONTHJan]"
## [65] "r_BIRD[nest-box2,MONTHJan]" "r_BIRD[nest-box3,MONTHJan]"
## [67] "r_BIRD[nest-box4,MONTHJan]" "r_BIRD[nest-box5,MONTHJan]"
## [69] "r_BIRD[nest-box6,MONTHJan]" "r_BIRD[nest-box7,MONTHJan]"
## [71] "r_BIRD[nest-box8,MONTHJan]" "r_BIRD[nest-box9,MONTHJan]"
## [73] "r_BIRD[other1,MONTHJan]" "r_BIRD[other10,MONTHJan]"
## [75] "r_BIRD[other2,MONTHJan]" "r_BIRD[other3,MONTHJan]"
## [77] "r_BIRD[other4,MONTHJan]" "r_BIRD[other5,MONTHJan]"
## [79] "r_BIRD[other6,MONTHJan]" "r_BIRD[other7,MONTHJan]"
## [81] "r_BIRD[other8,MONTHJan]" "r_BIRD[other9,MONTHJan]"
## [83] "r_BIRD[tree1,MONTHJan]" "r_BIRD[tree10,MONTHJan]"
## [85] "r_BIRD[tree2,MONTHJan]" "r_BIRD[tree3,MONTHJan]"
## [87] "r_BIRD[tree4,MONTHJan]" "r_BIRD[tree5,MONTHJan]"
## [89] "r_BIRD[tree6,MONTHJan]" "r_BIRD[tree7,MONTHJan]"
## [91] "r_BIRD[tree8,MONTHJan]" "r_BIRD[tree9,MONTHJan]"
## [93] "prior_Intercept" "prior_b"
## [95] "prior_sigma" "prior_sd_BIRD"
## [97] "prior_cor_BIRD" "lp__"
## [99] "accept_stat__" "stepsize__"
## [101] "treedepth__" "n_leapfrog__"
## [103] "divergent__" "energy__"
pars <- starling.brm3 %>% get_variables()
pars <- str_extract(pars, "^b_.*|^sigma$|^sd.*") %>% na.omit()
starling.brm3$fit %>%
stan_trace(pars = pars)
The chains appear well mixed and very similar
starling.brm3$fit %>%
stan_ac(pars = pars)
There is no evidence of auto-correlation in the MCMC samples
starling.brm3$fit %>% stan_rhat()
All Rhat values are below 1.05, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
starling.brm3$fit %>% stan_ess()
Ratios all very high.
starling.brm3$fit %>%
stan_dens(separate_chains = TRUE, pars = pars)
The ggmean package also as a set of MCMC diagnostic
functions. Lets start with the MCMC diagnostics.
Of these, we will focus on:
## starling.ggs <- starling.brm3 %>% ggs(burnin = FALSE, inc_warmup = FALSE)
## starling.ggs %>% ggs_traceplot()
The chains appear well mixed and very similar
## ggs_autocorrelation(starling.ggs)
There is no evidence of auto-correlation in the MCMC samples
## ggs_Rhat(starling.ggs)
All Rhat values are below 1.05, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
## ggs_effective(starling.ggs)
Ratios all very high.
## ggs_crosscorrelation(starling.ggs)
## ggs_grb(starling.ggs)
Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.
available_ppc()
## bayesplot PPC module:
## ppc_bars
## ppc_bars_grouped
## ppc_boxplot
## ppc_data
## ppc_dens
## ppc_dens_overlay
## ppc_dens_overlay_grouped
## ppc_ecdf_overlay
## ppc_ecdf_overlay_grouped
## ppc_error_binned
## ppc_error_hist
## ppc_error_hist_grouped
## ppc_error_scatter
## ppc_error_scatter_avg
## ppc_error_scatter_avg_vs_x
## ppc_freqpoly
## ppc_freqpoly_grouped
## ppc_hist
## ppc_intervals
## ppc_intervals_data
## ppc_intervals_grouped
## ppc_km_overlay
## ppc_loo_intervals
## ppc_loo_pit
## ppc_loo_pit_data
## ppc_loo_pit_overlay
## ppc_loo_pit_qq
## ppc_loo_ribbon
## ppc_ribbon
## ppc_ribbon_data
## ppc_ribbon_grouped
## ppc_rootogram
## ppc_scatter
## ppc_scatter_avg
## ppc_scatter_avg_grouped
## ppc_stat
## ppc_stat_2d
## ppc_stat_freqpoly_grouped
## ppc_stat_grouped
## ppc_violin_grouped
starling.brm3 %>% pp_check(type = "dens_overlay", nsamples = 100)
The model draws appear to be consistent with the observed data.
starling.brm3 %>% pp_check(type = "error_scatter_avg")
This is not really interpretable
starling.brm3 %>% pp_check(group = "BIRD", type = "intervals")
The shinystan package allows the full suite of MCMC
diagnostics and posterior predictive checks to be accessed via a web
interface.
# library(shinystan)
# launch_shinystan(starling.brm2)
DHARMa residuals provide very useful diagnostics. Unfortunately, we
cannot directly use the simulateResiduals() function to
generate the simulated residuals. However, if we are willing to
calculate some of the components yourself, we can still obtain the
simulated residuals from the fitted stan model.
We need to supply:
preds <- starling.brm3 %>% posterior_predict(nsamples = 250, summary = FALSE)
starling.resids <- createDHARMa(
simulatedResponse = t(preds),
observedResponse = starling$MASS,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = FALSE
)
plot(starling.resids, quantreg = TRUE)
Conclusions:
starling.brm3 %>%
conditional_effects("SITUATION:MONTH") %>%
plot(points = TRUE)
starling.brm3 %>%
ggpredict(~ SITUATION * MONTH) %>%
plot(add.data = TRUE)
starling.brm3 %>%
ggemmeans(~ SITUATION | MONTH) %>%
plot(add.data = TRUE)
Partial.obs <- starling.brm3$data %>%
mutate(
Pred = predict(starling.brm3)[, "Estimate"],
Resid = resid(starling.brm3)[, "Estimate"],
Obs = Pred + Resid
)
starling.brm3 %>%
fitted_draws(newdata = starling, re_formula = NA) %>%
median_hdci() %>%
ggplot(aes(x = SITUATION, y = .value, color = MONTH)) +
geom_pointrange(aes(ymin = .lower, ymax = .upper)) +
geom_line() +
geom_point(
data = Partial.obs, aes(y = Obs, x = SITUATION, color = MONTH),
position = position_nudge(x = 0.1)
) +
geom_point(
data = starling, aes(y = MASS, x = SITUATION, color = MONTH), alpha = 0.2,
position = position_nudge(x = 0.05)
)
brms captures the MCMC samples from stan
within the returned list. There are numerous ways to retrieve and
summarise these samples. The first three provide convenient numeric
summaries from which you can draw conclusions, the last four provide
ways of obtaining the full posteriors.
The summary() method generates simple summaries (mean,
standard deviation as well as 10, 50 and 90 percentiles).
starling.brm3 %>% summary()
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: MASS ~ MONTH * SITUATION + (MONTH | BIRD)
## Data: starling (Number of observations: 80)
## Draws: 3 chains, each with iter = 1000; warmup = 200; thin = 5;
## total post-warmup draws = 480
##
## Group-Level Effects:
## ~BIRD (Number of levels: 40)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.74 0.60 0.02 2.23 1.00 2061
## sd(MONTHJan) 1.34 1.08 0.05 3.84 1.00 1869
## cor(Intercept,MONTHJan) 0.03 0.58 -0.94 0.96 1.00 2314
## Tail_ESS
## sd(Intercept) 1992
## sd(MONTHJan) 1968
## cor(Intercept,MONTHJan) 2236
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 79.13 1.14 77.01 81.39 1.00 2354
## MONTHJan 8.53 1.53 5.48 11.53 1.00 2545
## SITUATIONnestMbox 0.38 1.59 -2.75 3.42 1.00 2577
## SITUATIONother -3.48 1.59 -6.53 -0.39 1.00 2489
## SITUATIONtree 4.15 1.60 1.06 7.27 1.00 2399
## MONTHJan:SITUATIONnestMbox 2.02 2.19 -2.24 6.44 1.00 2546
## MONTHJan:SITUATIONother 0.05 2.19 -4.24 4.49 1.00 2389
## MONTHJan:SITUATIONtree -0.95 2.17 -5.25 3.21 1.00 2474
## Tail_ESS
## Intercept 2335
## MONTHJan 2327
## SITUATIONnestMbox 2411
## SITUATIONother 2329
## SITUATIONtree 2259
## MONTHJan:SITUATIONnestMbox 2173
## MONTHJan:SITUATIONother 2370
## MONTHJan:SITUATIONtree 2428
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 3.97 0.41 3.17 4.81 1.00 2226 2145
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Conclusions:
starling.brm3$fit %>%
tidyMCMC(
estimate.method = "median",
conf.int = TRUE, conf.method = "HPDinterval",
rhat = TRUE, ess = TRUE
)
Conclusions:
see above
Due to the presence of a log transform in the predictor, it is better to use the regex version.
starling.brm3 %>% get_variables()
## [1] "b_Intercept" "b_MONTHJan"
## [3] "b_SITUATIONnestMbox" "b_SITUATIONother"
## [5] "b_SITUATIONtree" "b_MONTHJan:SITUATIONnestMbox"
## [7] "b_MONTHJan:SITUATIONother" "b_MONTHJan:SITUATIONtree"
## [9] "sd_BIRD__Intercept" "sd_BIRD__MONTHJan"
## [11] "cor_BIRD__Intercept__MONTHJan" "sigma"
## [13] "r_BIRD[inside1,Intercept]" "r_BIRD[inside10,Intercept]"
## [15] "r_BIRD[inside2,Intercept]" "r_BIRD[inside3,Intercept]"
## [17] "r_BIRD[inside4,Intercept]" "r_BIRD[inside5,Intercept]"
## [19] "r_BIRD[inside6,Intercept]" "r_BIRD[inside7,Intercept]"
## [21] "r_BIRD[inside8,Intercept]" "r_BIRD[inside9,Intercept]"
## [23] "r_BIRD[nest-box1,Intercept]" "r_BIRD[nest-box10,Intercept]"
## [25] "r_BIRD[nest-box2,Intercept]" "r_BIRD[nest-box3,Intercept]"
## [27] "r_BIRD[nest-box4,Intercept]" "r_BIRD[nest-box5,Intercept]"
## [29] "r_BIRD[nest-box6,Intercept]" "r_BIRD[nest-box7,Intercept]"
## [31] "r_BIRD[nest-box8,Intercept]" "r_BIRD[nest-box9,Intercept]"
## [33] "r_BIRD[other1,Intercept]" "r_BIRD[other10,Intercept]"
## [35] "r_BIRD[other2,Intercept]" "r_BIRD[other3,Intercept]"
## [37] "r_BIRD[other4,Intercept]" "r_BIRD[other5,Intercept]"
## [39] "r_BIRD[other6,Intercept]" "r_BIRD[other7,Intercept]"
## [41] "r_BIRD[other8,Intercept]" "r_BIRD[other9,Intercept]"
## [43] "r_BIRD[tree1,Intercept]" "r_BIRD[tree10,Intercept]"
## [45] "r_BIRD[tree2,Intercept]" "r_BIRD[tree3,Intercept]"
## [47] "r_BIRD[tree4,Intercept]" "r_BIRD[tree5,Intercept]"
## [49] "r_BIRD[tree6,Intercept]" "r_BIRD[tree7,Intercept]"
## [51] "r_BIRD[tree8,Intercept]" "r_BIRD[tree9,Intercept]"
## [53] "r_BIRD[inside1,MONTHJan]" "r_BIRD[inside10,MONTHJan]"
## [55] "r_BIRD[inside2,MONTHJan]" "r_BIRD[inside3,MONTHJan]"
## [57] "r_BIRD[inside4,MONTHJan]" "r_BIRD[inside5,MONTHJan]"
## [59] "r_BIRD[inside6,MONTHJan]" "r_BIRD[inside7,MONTHJan]"
## [61] "r_BIRD[inside8,MONTHJan]" "r_BIRD[inside9,MONTHJan]"
## [63] "r_BIRD[nest-box1,MONTHJan]" "r_BIRD[nest-box10,MONTHJan]"
## [65] "r_BIRD[nest-box2,MONTHJan]" "r_BIRD[nest-box3,MONTHJan]"
## [67] "r_BIRD[nest-box4,MONTHJan]" "r_BIRD[nest-box5,MONTHJan]"
## [69] "r_BIRD[nest-box6,MONTHJan]" "r_BIRD[nest-box7,MONTHJan]"
## [71] "r_BIRD[nest-box8,MONTHJan]" "r_BIRD[nest-box9,MONTHJan]"
## [73] "r_BIRD[other1,MONTHJan]" "r_BIRD[other10,MONTHJan]"
## [75] "r_BIRD[other2,MONTHJan]" "r_BIRD[other3,MONTHJan]"
## [77] "r_BIRD[other4,MONTHJan]" "r_BIRD[other5,MONTHJan]"
## [79] "r_BIRD[other6,MONTHJan]" "r_BIRD[other7,MONTHJan]"
## [81] "r_BIRD[other8,MONTHJan]" "r_BIRD[other9,MONTHJan]"
## [83] "r_BIRD[tree1,MONTHJan]" "r_BIRD[tree10,MONTHJan]"
## [85] "r_BIRD[tree2,MONTHJan]" "r_BIRD[tree3,MONTHJan]"
## [87] "r_BIRD[tree4,MONTHJan]" "r_BIRD[tree5,MONTHJan]"
## [89] "r_BIRD[tree6,MONTHJan]" "r_BIRD[tree7,MONTHJan]"
## [91] "r_BIRD[tree8,MONTHJan]" "r_BIRD[tree9,MONTHJan]"
## [93] "prior_Intercept" "prior_b"
## [95] "prior_sigma" "prior_sd_BIRD"
## [97] "prior_cor_BIRD" "lp__"
## [99] "accept_stat__" "stepsize__"
## [101] "treedepth__" "n_leapfrog__"
## [103] "divergent__" "energy__"
starling.draw <- starling.brm3 %>%
gather_draws(`b.Intercept.*|b_SITUATION.*|b_MONTH.*`, regex = TRUE)
starling.draw
We can then summarise this
starling.draw %>% median_hdci()
starling.brm3 %>%
gather_draws(`b_Intercept.*|b_SITUATION.*|b_MONTH.*`, regex = TRUE) %>%
ggplot() +
geom_vline(xintercept = 0, linetype = "dashed") +
stat_slab(aes(
x = .value, y = .variable,
fill = stat(ggdist::cut_cdf_qi(cdf,
.width = c(0.5, 0.8, 0.95),
labels = scales::percent_format()
))
), color = "black") +
scale_fill_brewer("Interval", direction = -1, na.translate = FALSE)
starling.brm3 %>%
gather_draws(`.Intercept.*|b_SITUATION.*|b_MONTH.*`, regex = TRUE) %>%
ggplot() +
geom_vline(xintercept = 0, linetype = "dashed") +
stat_halfeye(aes(x = .value, y = .variable)) +
theme_classic()
starling.brm3$fit %>% plot(type = "intervals")
This is purely a graphical depiction on the posteriors.
starling.brm3 %>% tidy_draws()
starling.brm3 %>% spread_draws(`.*Intercept.*|b_SITUATION.*|b_MONTH.*`, regex = TRUE)
starling.brm3 %>%
posterior_samples() %>%
as_tibble()
starling.brm3 %>%
bayes_R2(re.form = NA, summary = FALSE) %>%
median_hdci()
starling.brm3 %>%
bayes_R2(re.form = ~ (1 | BIRD), summary = FALSE) %>%
median_hdci()
starling.brm3 %>%
bayes_R2(re.form = ~ (MONTH | BIRD), summary = FALSE) %>%
median_hdci()
starling.brm3 %>%
emmeans(~SITUATION) %>%
pairs()
## contrast estimate lower.HPD upper.HPD
## inside - (nest-box) -1.41 -3.783 1.182
## inside - other 3.48 0.759 5.881
## inside - tree -3.68 -6.203 -1.037
## (nest-box) - other 4.86 2.321 7.406
## (nest-box) - tree -2.30 -5.026 0.263
## other - tree -7.14 -9.789 -4.419
##
## Point estimate displayed: median
## HPD interval probability: 0.95
## contrast estimate lower.HPD upper.HPD
## Nov - Jan -8.81 -10.6 -7.08
##
## Point estimate displayed: median
## HPD interval probability: 0.95
## SITUATION = inside:
## contrast estimate lower.HPD upper.HPD
## Nov - Jan -8.56 -11.4 -5.38
##
## SITUATION = nest-box:
## contrast estimate lower.HPD upper.HPD
## Nov - Jan -10.54 -14.1 -7.22
##
## SITUATION = other:
## contrast estimate lower.HPD upper.HPD
## Nov - Jan -8.55 -11.8 -4.91
##
## SITUATION = tree:
## contrast estimate lower.HPD upper.HPD
## Nov - Jan -7.58 -11.0 -4.29
##
## Point estimate displayed: median
## HPD interval probability: 0.95
Compare specific sets of SITUATION within each MONTH
levels(starling$SITUATION)
## [1] "inside" "nest-box" "other" "tree"
## MONTH = Nov:
## contrast estimate lower.HPD upper.HPD
## SITUATION.Comp1 3.666 0.855 6.49
## SITUATION.Comp2 1.544 -1.131 4.11
##
## MONTH = Jan:
## contrast estimate lower.HPD upper.HPD
## SITUATION.Comp1 4.610 1.401 7.69
## SITUATION.Comp2 0.552 -2.590 3.46
##
## Point estimate displayed: median
## HPD interval probability: 0.95